########################
# Plot Counterfactuals #
########################
pdf(paste("DFQF",str.sample,str.rho,str.cb,str.robust,str.boot,".pdf",sep=""));

#par(par.current);
# DF #
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Latent and Counterfactual Distributions, ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.F1,rev(lbound.F1)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(ys, F1, lwd=1, col = 4  );


polygon(c(ys,rev(ys)),c(ubound.F2_c,rev(lbound.F2_c)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(ys, F2_c, lwd=1, col = 'dark green'  );

polygon(c(ys,rev(ys)),c(ubound.F2,rev(lbound.F2)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(ys, F2, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter1,paste("Fz=",legend.counter2,",Theta=",legend.counter1,sep=""),legend.counter2),
       col= c(4,"dark green","dark grey"),lty=1);

# DF with Joint CB #
plot(range(ys), c(0,1), type="n",col=1, xlab="Log of hourly wage", lwd=2, ylab="Probability", 
     main=paste("Latent and Counterfactual Distributions with Joint CB, ",main.str, sep=""), sub=" ");

polygon(c(ys,rev(ys)),c(ubound.F.joint[,1],rev(lbound.F.joint[,1])), density=100, border=F, col='light blue', lty=1, lwd=1);
lines(ys, F1, lwd=1, col = 4  );

polygon(c(ys,rev(ys)),c(ubound.F.joint[,2],rev(lbound.F.joint[,2])), density=100, border=F, col='light green', lty=1, lwd=1);
lines(ys, F2_c, lwd=1, col = 'dark green'  );

polygon(c(ys,rev(ys)),c(ubound.F.joint[,3],rev(lbound.F.joint[,3])), density=100, border=F, col='light grey', lty=1, lwd=1);
lines(ys, F2, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c(legend.counter1,paste("Fz=",legend.counter2,",Theta=",legend.counter1,sep=""),legend.counter2),
       col= c(4,"dark green","dark grey"),lty=1);


# QF #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Latent and Counterfactual Quantile Functions, ",main.str, sep=""), sub=" ");

polygon(c(ubound.F1,rev(lbound.F1)),c(ys,rev(ys)), density=100, border=F, col='light blue', lty=1, lwd=1);
lines(F1, ys, lwd=1, col = 4  );

polygon(c(ubound.F2_c,rev(lbound.F2_c)),c(ys,rev(ys)), density=100, border=F, col='light green', lty=1, lwd=1);
lines(F2_c, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.F2,rev(lbound.F2)),c(ys,rev(ys)), density=100, border=F, col='light grey', lty=1, lwd=1);
lines(F2, ys, lwd=1, col = 'dark grey'  );


legend("bottomright",legend=c(legend.counter1,paste("Fz=",legend.counter2,",Theta=",legend.counter1,sep=""),legend.counter2),
       col= c(4,"dark green","dark grey"),lty=1);


# QF with Joint CB #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Latent and Counterfactual Quantile Functions with Joint CB, ",main.str, sep=""), sub=" ");

polygon(c(ubound.F.joint[,1],rev(lbound.F.joint[,1])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F1, ys, lwd=1, col = 4  );

polygon(c(ubound.F.joint[,2],rev(lbound.F.joint[,2])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2_c, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.F.joint[,3],rev(lbound.F.joint[,3])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F2, ys, lwd=1, col = 'dark grey'  );


legend("bottomright",legend=c(legend.counter1,paste("Fz=",legend.counter2,",Theta=",legend.counter1,sep=""),legend.counter2),
       col= c(4,"dark green","dark grey"),lty=1);


# QF with Bootstrap CB #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2,
     ylab="Log of hourly wage", 
     main=paste("Latent and Counterfactual Quantile Functions with Joint Bootstrap CB, ",main.str, sep=""),
     sub=" ");

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),c(ubound.QF[,1],rev(lbound.QF[,1])), density=100, border=F,
        col='light blue', lty = 1, lwd = 1);
lines(F1, ys, lwd=1, col = 4  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),c(ubound.QF[,2],rev(lbound.QF[,2])), density=100, border=F,
        col='light green', lty = 1, lwd = 1);
lines(F2_c, ys, lwd=1, col = 'dark green'  );

polygon(c(seq(l.thres,u.thres,by.thres),seq(u.thres,l.thres,-by.thres)),c(ubound.QF[,3],rev(lbound.QF[,3])), density=100, border=F,
        col='light grey', lty = 1, lwd = 1);
lines(F2, ys, lwd=1, col = 'dark grey'  );


legend("bottomright",legend=c(legend.counter1,paste("Fz=",legend.counter2,",Theta=",legend.counter1,sep=""),legend.counter2),
       col= c(4,"dark green","dark grey"),lty=1);


# QF with Average Effects #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Latent and Counterfactual Quantile Functions, ",main.str, sep=""), sub=" ");

polygon(c(ubound.F1,rev(lbound.F1)),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F1, ys, lwd=1, col = 4  );

polygon(c(ubound.F2_c,rev(lbound.F2_c)),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2_c, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.F2,rev(lbound.F2)),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F2, ys, lwd=1, col = 'dark grey'  );

abline(h=e.avg.F1, col=4);
abline(h=ubound.avg[1], col=4, lty=2);
abline(h=lbound.avg[1], col=4, lty=2);

abline(h=e.avg.F2_c, col='dark green');
abline(h=ubound.avg[2], col='dark green', lty=2);
abline(h=lbound.avg[2], col='dark green', lty=2);

abline(h=e.avg.F2, col='dark grey');
abline(h=ubound.avg[3], col='dark grey', lty=2);
abline(h=lbound.avg[3], col='dark grey', lty=2);


legend("bottomright",legend=c(legend.counter1,paste("Fz=",legend.counter2,",Theta=",legend.counter1,sep=""),legend.counter2),
       col= c(4,"dark green","dark grey"),lty=1);

# QF with Average Effects and joint CB #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Latent and Counterfactual Quantile Functions with Joint CB, ",main.str, sep=""), sub=" ");

polygon(c(ubound.F.joint[,1],rev(lbound.F.joint[,1])),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F1, ys, lwd=1, col = 4  );

polygon(c(ubound.F.joint[,2],rev(lbound.F.joint[,2])),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2_c, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.F.joint[,3],rev(lbound.F.joint[,3])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F2, ys, lwd=1, col = 'dark grey'  );

abline(h=e.avg.F1, col=4);
abline(h=ubound.avg.joint[1], col=4, lty=2);
abline(h=lbound.avg.joint[1], col=4, lty=2);

abline(h=e.avg.F2_c, col='dark green');
abline(h=ubound.avg.joint[2], col='dark green', lty=2);
abline(h=lbound.avg.joint[2], col='dark green', lty=2);

abline(h=e.avg.F2, col='dark grey');
abline(h=ubound.avg.joint[3], col='dark grey', lty=2);
abline(h=lbound.avg.joint[3], col='dark grey', lty=2);


legend("bottomright",legend=c(legend.counter1,paste("Fz=",legend.counter2,",Theta=",legend.counter1,sep=""),legend.counter2),
       col= c(4,"dark green","dark grey"),lty=1);


# QF with Heckman Decomposition #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Heckman Latent and Counterfactual Quantile Functions, ",main.str, sep=""), sub=" ");

polygon(c(ubound.F1,rev(lbound.F1)),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F1, ys, lwd=1, col = 4  );

polygon(c(ubound.F2_c,rev(lbound.F2_c)),c(ys,rev(ys)), density=100, border=F, col='light green', lty = 1, lwd = 1);
lines(F2_c, ys, lwd=1, col = 'dark green'  );

polygon(c(ubound.F2,rev(lbound.F2)),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F2, ys, lwd=1, col = 'dark grey'  );

lines(F.Heck1, ys, lwd=1, col=4);
lines(F.Heck2_c, ys, lwd=1, col='dark green');
lines(F.Heck2, ys, lwd=1, col='dark grey');

legend("bottomright",legend=c(legend.counter1,paste("Fz=",legend.counter2,",Theta=",legend.counter1,sep=""),legend.counter2),
       col= c(4,"dark green","dark grey"),lty=1);


# Percentage Change #
curve(pct.F.theta, ylim=ylim.pct.F, col=2, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage Decomposition of Distribution Change", from=l.trim/100,to=u.trim/100);
curve(lbound.pct.F.theta, add=T, col=2, lty=2);
curve(ubound.pct.F.theta, add=T, col=2, lty=2);
curve(pct.Heck.F.theta, add=T, col=1, lty=2);
legend("bottomleft", legend=c("QF %Decomposition","95% CB","Heckman"), col=c(2,2,1),lty=c(1,2,2));


# Percentage Change of DF #
matplot(ys[qn.plot], cbind(pct.DF.theta, ubound.pct.DF.theta, lbound.pct.DF.theta, pct.Heck.DF.theta)[qn.plot,], 
        col=c(2,2,2,1), type="l", lty=c(1,2,2,2),
        ylab="% of Decomposition", xlab="Quantile", 
        main="Percentage Decomposition of Distribution Change -- DF on 5%-95%");
legend("bottomleft", legend=c("DF %Decomposition","95% CB","Heckman"), col=c(2,2,1),lty=c(1,2,2));


# Quantile Function Change #
curve(QF2_c(x)-QF2(x), ylim=ylim.QF.diff, ylab="Decomposition Component", xlab="Quantile", 
      main="Decomposition of Distribution Change", col=2, from=l.trim/100, to=u.trim/100);
curve(ubound.QF.joint$F2_c(x)-lbound.QF.joint$F2(x), add=T, col=2, lty=2);
curve(lbound.QF.joint$F2_c(x)-ubound.QF.joint$F2(x), add=T, col=2, lty=2);
curve(QF.Heck2_c(x) - QF.Heck2(x), add=T, col="dark red", lty=3);
curve(QF1(x)-QF2_c(x), add=T, col=3);
curve(ubound.QF.joint$F1(x)-lbound.QF.joint$F2_c(x), add=T, col=3, lty=2);
curve(lbound.QF.joint$F1(x)-ubound.QF.joint$F2_c(x), add=T, col=3, lty=2);
curve(QF.Heck1(x) - QF.Heck2_c(x), add=T, col="dark green", lty=3);
curve(QF1(x)-QF2(x), add=T, col=1);
curve(ubound.QF.joint$F1(x)-lbound.QF.joint$F2(x), add=T, col=1, lty=2);
curve(lbound.QF.joint$F1(x)-ubound.QF.joint$F2(x), add=T, col=1, lty=2);
curve(QF.Heck1(x) - QF.Heck2(x), add=T, col=1, lty=3);
abline(h=0);
legend("bottomright", legend=c(as.expression(bquote("QF2|"*beta*"=1 - QF2")),
                               as.expression(bquote("QF1 - QF2|"*beta*"=1")), 
                               as.expression(bquote("QF1 - QF2"))), col=c(2,3,1),lty=1);

# Percentage Change with Joint CB #
curve(pct.F.theta, ylim=ylim.pct.F, col=2, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage Decomposition of Distribution Change with Joint CB", from=l.trim/100,to=u.trim/100);
curve((ubound.QF.joint$F2_c(x)-lbound.QF.joint$F2(x))/(lbound.QF.joint$F1(x)-ubound.QF.joint$F2(x))*100, add=T, col=2, lty=2);
curve((lbound.QF.joint$F2_c(x)-ubound.QF.joint$F2(x))/(ubound.QF.joint$F1(x)-lbound.QF.joint$F2(x))*100, add=T, col=2, lty=2);
curve(pct.Heck.F.theta, add=T, col=1, lty=2);
legend("bottomleft", legend=c("QF %Decomposition","95% CB","Heckman"), col=c(2,2,1),lty=c(1,2,2));


# Comparison of Distributions #
# Counter Level 1 -- F vs. F.cond #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("QF in Full and Observed Populations, ",main.str.counter1, sep=""), sub=" ");

polygon(c(ubound.F1,rev(lbound.F1)),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F1, ys, lwd=1, col = 4  );

polygon(c(ubound.cond[,5],rev(lbound.cond[,5])), c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col='dark grey');

legend("bottomright",legend=c("Latent","Observed"),col= c(4,'dark grey'),lty=1);


# Counter Level 1 -- F.cond vs. F.obs #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("QF in Observed Population and Empirical QF, ",main.str.counter1, sep=""), sub=" ");
polygon(c(ubound.cond[,5],rev(lbound.cond[,5])), c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1.cond, ys, lwd=1, col='dark grey');
lines(F1.obs, ys, lwd=1, col=1);

legend("bottomright",legend=c("Observed","Empirical"), col= c('dark grey',1),lty=1);


# Counter Level 2 -- F vs. F.cond #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("QF in Full and Observed populations, ",main.str.counter2, sep=""), sub=" ");

polygon(c(ubound.F2,rev(lbound.F2)),c(ys,rev(ys)), density=100, border=F, col='light blue', lty = 1, lwd = 1);
lines(F2, ys, lwd=1, col = 4  );

polygon(c(ubound.cond[,1],rev(lbound.cond[,1])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'dark grey'  );

legend("bottomright",legend=c("Latent","Observed"), col= c(4, 'dark grey'),lty=1);


# Counter Level 2 -- F.cond vs. F.obs #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("QF in Observed Population and Empirical QF, ",main.str.counter2, sep=""), sub=" ");

polygon(c(ubound.cond[,1],rev(lbound.cond[,1])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F2.cond, ys, lwd=1, col = 'dark grey'  );
lines(F2.obs, ys, lwd=1, col=1);
legend("bottomright",legend=c("Observed","Empirical"), col= c('dark grey',1),lty=1);



##################
# Plot for paper #
##################
# QF #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main="Offered wage");

polygon(c(ubound.F1,rev(lbound.F1)),c(ys,rev(ys)), density=100, border=F, col='light grey', lty=1, lwd=1);
lines(F1, ys, lwd=1, col = 'dark grey' );

polygon(c(ubound.F2,rev(lbound.F2)),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty=1, lwd=1);
lines(F2, ys, lwd=1, col = 'firebrick'  );

legend("bottomright",legend=c(legend.counter2,legend.counter1), col= c("firebrick","dark grey"),lty=1);


# QF with Joint CB #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Latent QF with Joint CB, ",main.str, sep=""), sub=" ");

polygon(c(ubound.F.joint[,1],rev(lbound.F.joint[,1])),c(ys,rev(ys)), density=100, border=F, col='light grey', lty = 1, lwd = 1);
lines(F1, ys, lwd=1, col = 'dark grey'  );

polygon(c(ubound.F.joint[,3],rev(lbound.F.joint[,3])),c(ys,rev(ys)), density=100, border=F, col='firebrick1', lty = 1, lwd = 1);
lines(F2, ys, lwd=1, col = 'firebrick'  );

legend("bottomright",legend=c(legend.counter2,legend.counter1), col= c("firebrick","dark grey"),lty=1);


# Percentage Change #
curve(pct.F.theta, ylim=ylim.pct.F, col=4, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage Decomposition of Latent QF", from=l.trim/100,to=u.trim/100);
curve(lbound.pct.F.theta, add=T, col=4, lty=2);
curve(ubound.pct.F.theta, add=T, col=4, lty=2);
curve(pct.Heck.F.theta, add=T, col="dark blue", lty=3);
abline(h=0);
legend("bottomright", legend=c(expression(beta)), col=4, lty=1);


# Percentage Change with Joint CB #
curve(pct.F.theta, ylim=ylim.pct.F, col=4, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage Decomposition of Latent QF with Joint CB", from=l.trim/100,to=u.trim/100);
curve((ubound.QF.joint$F2_c(x)-lbound.QF.joint$F2(x))/(lbound.QF.joint$F1(x)-ubound.QF.joint$F2(x))*100, add=T, col=4, lty=2);
curve((lbound.QF.joint$F2_c(x)-ubound.QF.joint$F2(x))/(ubound.QF.joint$F1(x)-lbound.QF.joint$F2(x))*100, add=T, col=4, lty=2);
curve(pct.Heck.F.theta, add=T, col='dark blue', lty=3);
abline(h=0);
legend("bottomright", legend=c(expression(beta)), col=4,lty=1);


# Percentage Change -- Decomposition by two parts #
curve(pct.F.theta, ylim=ylim.pct.F, col=4, ylab="% of Decomposition", xlab="Quantile Index",
      main="Percentage decomposition", from=l.trim/100,to=u.trim/100);
curve(lbound.pct.F.theta, add=T, col=4, lty=2);
curve(ubound.pct.F.theta, add=T, col=4, lty=2);
curve(pct.Heck.F.theta, add=T, col="dark blue", lty=3);
curve(100 - pct.F.theta(x), add=T, col=6)
curve(100 - lbound.pct.F.theta(x), add=T, col=6, lty=2);
curve(100 - ubound.pct.F.theta(x), add=T, col=6, lty=2);
curve(100 - pct.Heck.F.theta(x), add=T, col=6, lty=3);
abline(h=0);
legend("bottomright", legend=c(expression(beta), "Fx"), col=c(4,6), lty=1);


dev.off()


pdf(paste("Comp_Dist", str.sample,str.rho,str.cb,str.robust,str.boot,".pdf", sep = ""))

# Comparison of Distributions #
# Counter Level 1 -- F vs. F.cond #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=main.str.counter1);

lines(F1, ys, lwd=1, col = 4  );

lines(F1.cond, ys, lwd=1, col='orange');

legend("bottomright",legend=c("Latent","Observed"),col= c(4,'orange'),lty=1);

# Counter Level 2 -- F vs. F.cond #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=main.str.counter2);

lines(F2, ys, lwd=1, col = 4  );

lines(F2.cond, ys, lwd=1, col='orange');

legend("bottomright",legend=c("Latent","Observed"), col= c(4, 'orange'),lty=1);


# Counter Level 1 -- F.cond vs. F.obs #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Model and Empirical QF in Observed Population, ",main.str.counter1, sep=""), sub=" ");
lines(F1.cond, ys, lwd=3, col='dark grey');
lines(F1.obs, ys, lwd=1, col=1);

legend("bottomright",legend=c("Model","Empirical"), col= c('dark grey',1),lty=1);


# Counter Level 2 -- F.cond vs. F.obs #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("Model and Empirical QF in Observed Population, ",main.str.counter2, sep=""), sub=" ");
lines(F2.cond, ys, lwd=3, col = 'dark grey'  );
lines(F2.obs, ys, lwd=1, col=1);
legend("bottomright",legend=c("Model","Empirical"), col= c('dark grey',1),lty=1);

# Counter Level 1 -- F.cond vs. F.noselection #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("QF in Observed Population with and without Selection, ",main.str.counter1, sep=""), sub=" ");
lines(F1.noselection, ys, lwd=1, col=4);
lines(F1.cond, ys, lwd=1, col='orange');

legend("bottomright",legend=c("Without Selection","With Selection"), col= c(4,'orange'),lty=1);


# Counter Level 2 -- F.cond vs. F.noselection #
plot(c(0,1), range(ys), type="n",col=1, xlab="Quantile Index", lwd=2, ylab="Log of hourly wage", 
     main=paste("QF in Observed Population with and without Selection, ",main.str.counter2, sep=""), sub=" ");
lines(F2.noselection, ys, lwd=1, col = 4  );
lines(F2.cond, ys, lwd=1, col='orange');
legend("bottomright",legend=c("Without Selection","With Selection"), col= c(4,'orange'),lty=1);


dev.off()
